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New  and  Improved  Solar  Radiation  Models  for  GPS  Satellites  Based  on  Flight 

Data 


Executive  summary 

This  report  documents  the  analysis  and  results  from  a  year-long  study  of  a  novel  approach  to  the 
problem  of  developing  solar  radiation  models  for  GPS  satellites.  The  approach  is  aiming  to 
replace  the  pre-launch  design  phase  of  solar  pressure  and  heat  reradiation  models  by  a  less  costly 
and  more  accurate  post-launch  phase.  The  approach  is  also  suitable  for  many  other  Earth-orbiting 
satellites.  In  this  approach  we  exploit  the  fact  that  a  collection  of  individual  orbit  solutions 
contains  more  information  about  the  dynamics  of  the  satellite  than  the  information  that  was  used 
to  generate  the  individual  solutions. 

The  current  GPS  constellation  of  Block  II  and  Block  IIA  satellites  was  used  as  a  prototype  for 
developing  and  validating  our  approach.  We  have  used  daily  GPS  precise  ephemerides 
(produced  routinely  at  JPL  for  the  International  GPS  Service  (IGS))  over  a  period  of  9  months  to 
adjust  a  parametrized  model  of  the  solar  pressure  so  as  to  obtain  best  fit.  The  resulting  model 
proved  to  be  more  accurate  than  the  standard  solar  pressure  model  for  GPS  satellites  (Fliegel  et 
al.,  1992).  In  a  separate  effort  we  developed  from  first  principles  a  new  solar  pressure  model  for 
eclipsing  GPS  satellites.  This  model  also  represents  a  significant  improvement  over  the  standard 
model. 

The  main  benefits  from  the  new  approach  are: 

-  Cost  reduction  due  to  the  replacement  of  costly  pre-launch  design  phase  with  cheap  post-launch 
process. 

-  Increased  accuracy  due  to  sensing  of  actual  satellite  behavior,  as  opposed  to  theoretical 
behavior. 

-  Increased  accuracy  by  allowing  an  infinite  fine-tuning  process. 

Increased  model  accuracy,  in  turn,  will  improve: 

-  Orbit  prediction  capability  including  autonav. 

-  Improved  UT 1  estimation  capabilities 

-  Improved  accuracy  for  any  high-precision  GPS  application,  in  particular,  scientific  applications 

such  as  geodesy,  meteorology  and  altimetry. 

-  Increase  our  understanding  of  the  space  environment  and  in-orbit  satellite  behavior. 

The  most  immediate  application  of  the  new  approach  is  for  the  refinement  of  the  nominal  solar 
radiation  model  for  Block  HR  satellites  (Fliegel  and  Gallini,  1997)  immediately  after  launch. 
Model  refinement  can  be  made  monthly  and  it  is  expected  that  each  refinement  will  be 
substantial  for  at  least  the  first  year,  given  the  typical  BOF  anomalies.  The  greatest  utility  will  be 
achieved  by  applying  the  new  approach  to  the  Block  IIF  satellites.  A  significant  saving  in  cost 
can  be  made  by  eliminating  the  pre-launch  design  phase  of  the  solar  radiation  model.  In  addition, 
lessons  learned  from  the  analysis  of  the  Block  II  solar  pressure  can  be  applied  to  the  design  of 
the  attitude  control  system  of  Block  IIF  satellites  and  improve  the  accuracy  of  any  solar  radiation 
model  developed  for  these  satellites. 

The  work  presented  here  is  a  proof  of  concept.  Although  a  significant  improvement  in  the  solar 
radiation  model  for  Block  II/IIA  satellites  was  achieved,  the  case  was  treated  as  a  prototype  and 
it  was  impossible  in  the  limited  time  we  had  to  explore  all  the  possible  modeling  options. 
Instead,  most  of  the  efforts  were  spent  overcoming  the  various  technical  difficulties  in 
implementing  our  approach  and  on  refining  and  tuning  the  technology.  Consequently,  further 
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accuracy  improvements  to  Block  II/IIA  solar  radiation  models  can  be  realized  with  only  a 
modest  amount  of  additional  analysis. 
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Introduction 

Orbiting  at  an  altitude  of  about  20,000  km,  with  no  drag  and  with  limited  sensitivity  to  the  details 
of  the  Earth’s  gravitational  pull,  the  GPS  satellites  seem  in  little  need  for  a  complex  dynamical 
model.  Yet  the  relatively  poor  geometry  of  the  observation  system  (providing  mostly  radial 
information),  combined  with  the  demand  of  some  applications  for  extremely  high  accuracy, 
create  the  need  for  a  very  careful  modeling  of  the  forces  acting  on  a  GPS  satellite. 

Of  those  forces,  solar  radiation  is  the  most  delicate  to  model  (Fliegel  et  al.,  1992).  Although  the 
radiation  force,  at  GPS  altitude,  is  only  second  in  magnitude  to  the  gravitational  pull  from  the 
Sun,  Earth  and  Moon,  its  uncertainties  are  much  higher  (Table  1).  The  solar  radiation  model  has 
received  relatively  little  attention  since  the  publication  of  the  T20  model  by  Fliegel  et  al.  (1992) 
and  is  considered  to  be  the  largest  error  source  in  GPS  orbit  determination  (Bertiger  et.  al.  1994, 
Gold  et.  al.  1995,  Beutler  et  al.,  1994).  The  errors  associated  with  the  solar  pressure  model  are 
particularly  large  during  the  satellite’s  eclipse  period  (Bar-Sever  et  al.,  1996,  Watkins  et  al., 
1996,  Fliegel  et  al.,  1992). 


Source 

Magnitude,  m 

Uncertainty,  % 

J2 

15,600 

0.00001 

Moon 

880 

small 

Sun  (gravitation) 

420 

small 

Solar  pressure 

130 

1-5 

C(2,2),  S(2,2) 

120 

0.001 

C(3,m).  S(3,m)  all  m 

20 

0.01 

Table  1.  The  perturbations  on  a  GPS  satellite  after  one  revolution  and  their  uncertainties.  The 
uncertainties  in  the  gravitational  forces  are  based  on  JGM3  formal  errors.  The  magnitude  of  the 
perturbations  are  based  on  Fliegel  and  Gallini  (1996). 


The  application  most  sensitive  to  force  modeling  errors  is  orbit  prediction.  This  is  because,  in  the 
absence  of  constraining  measurements,  the  orbit  will  typically  diverge  quadratically  from  the 
truth  (Colombo,  1989).  Certain  limits  on  GPS  orbit  prediction  errors  are  required  by  the  GPS 
operators  at  the  US  Air  Force  for  proper  maintenance  of  the  constellation.  Block  HR  satellites 
will  have  “autonav”  capabilities  which  require  quality  long-term  orbit  prediction  to  allow  for  the 
positioning  of  the  satellites  in  the  absence  of  up-link  communications  over  periods  of  several 
months.  Many  of  the  recently-emerged  real-time  applications  of  the  GPS,  such  as  WAAS  and 
weather  forecasting,  are  dependent  on  the  quality  of  GPS  orbit  prediction  (Bar-Sever,  1996a, 
Bertiger,  1997).  All  these  applications  stand  to  benefit  from  improvements  to  the  force  models. 

As  mentioned  above,  solar  radiation  mismodeling  is  a  limiting  error  source  in  precise 
applications  such  as  geodetic  positioning  and  low  earth  orbiter  tracking.  In  particular,  systematic 
errors  in  modeling  the  solar  radiation  alias  to  GPS -based  UT1  rate  estimates  through  the 
coupling  effect  of  the  right  ascension  of  the  ascending  node.  This  can  cause  a  uniform  rotation  of 
the  GPS  constellation  to  be  aliased  into  the  UT1  rate  estimate  (Fichten  et  al.,  1992).  Such  a 
uniform  rotation  can  be  caused  by  radiation  forces  such  as  the  Y  bias.  Reducing  the  uncertainties 
in  the  radiation  force  model,  especially  the  Y  bias,  will  contribute  to  reduce  this  aliasing  and 
improve  the  GPS -based  estimates  of  UT1  rate. 
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The  development  of  a  solar  radiation  model  for  a  spacecraft  is  usually  carried  out  by  the 
spacecraft  manufacturer  and  is  part  of  their  overall  contractual  obligation.  It  is  to  be  completed 
before  the  spacecraft  is  launched.  This  was  the  case  for  the  Block  I,  II  and  IIA  satellites  for 
which  the  ROCK4  and  ROCK42  solar  pressure  models  were  developed  by  Rockwell,  and  it  is 
also  the  case  for  the  Block  HR  satellites  developed  by  Lockheed  Martin.  The  cost  of  the 
development  of  a  solar  radiation  model  can  reach  millions  of  dollars,  depending  on  the  type  of 
the  mission.  In  the  case  of  the  Block  I,  II  and  IIA  satellites,  the  ROCK  models  (ROCK4  for 
Block  I  and  ROCK42  for  block  II  and  IIA)  were  developed  by  Rockwell,  programmed  by  IBM 
Federal  Systems  and  revised  and  improved  by  Rockwell  and  The  Aerospace  Corporation  (Fliegel 
and  Gallini,  1996).  The  best  model  currently  available  for  Block  II/IIA  satellites  is  the  T20 
model,  which  is  a  Fourier  series  representation  of  a  revised  ROCK42  model  (Fliegel  et  al., 
1992). 

The  ground-based  design  process  (as  opposed  to  the  flight-data-based  design  process  we  are 
proposing)  is  carried  out  by  modeling  the  spacecraft  as  a  collection  of  components,  each  with  its 
own  shape,  size,  optical  and  thermal  characteristics.  Given  the  nominal  mission  profile  and 
orbital  geometry,  the  process  employs  various  ray  tracing,  finite  element  and  finite  difference 
techniques  to  simulate  the  effects  of  impinging  photons  on  the  spacecraft,  and  derive  the  solar 
radiation  model.  This  is  an  exceedingly  complicated  process  because  most  spacecrafts  have  a 
complicated  shape  and  are  made  of  many  different  materials.  As  a  results,  it  is  necessary  to  make 
some  simplifying  assumptions  in  the  simulation.  For  example,  some  small  sub-structures  are 
ignored,  optical  and  thermal  properties  are  approximated  (rather  crudely),  mutual  shadowing  of 
components  is  applied  only  to  the  most  dominant  structures  (if  at  all),  and  secondary  reflections 
(from  one  component  to  the  other)  are  usually  ignored.  The  result  is  usually  a  fairly  good  model  - 
if  the  spacecraft's  in-orbit  behavior  does  not  deviate  much  from  the  nominal.  But  the  deficiencies 
of  such  an  approach  are  clear: 

1.  The  in-orbit  satellite  behavior  may  deviate  from  the  nominal.  Misorientation,  bending  and 
flexing  of  structures  are  quite  common.  For  example,  the  non-nominal  Y  bias  force  is  attributed 
to  solar  array  misalignment  (Fliegel  et  al.,  1992). 

2.  It  is  impossible  to  gauge  the  combined  effects  of  all  the  approximations,  simplifying 
assumptions  and  outright  errors  that  went  into  the  model.  Hence  the  actual  accuracy  of  the  model 
can  only  be  roughly  estimated. 

3.  The  model  is  not  adjustable  or  tunable.  If  accuracy  requirements  change  during  the  lifetime 
of  the  mission  it  usually  requires  a  costly  redesign  process. 

A  flight-data-based  design  process,  in  contrast,  has  the  following  advantages: 

1.  It  is  cheap,  since  all  the  data  is  available  and  is  in  conventional  form,  and  the  models  are 
simple  (Fourier  expansions). 

2.  It  reflects  actual  in-orbit  behavior. 

3.  It  is  more  accurate.  It  directly  accounts  for  the  combined  radiation  pressure  of  all  spacecraft 
components. 

4.  It  includes  radiation  effects  due  to  Earth  albedo. 

5.  It  is  infinitely  tunable  and  adjustable,  i.e.,  once  the  design  process  is  set  up  it  can  be  carried 
out  indefinitely  and  continuously  improve  and  adjust  for  changing  spacecraft  and  environmental 
conditions. 
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6.  It  provides  a  tool  for  learning  about  actual  in-orbit  behavior  of  satellites  and  for  flagging 
and  monitoring  problem  satellites. 

In  the  following  sections  we  will  describe  our  approach  to  the  design  process  of  flight-data-based 
solar  radiation  models  for  GPS  satellites,  and  we  will  present  evidence  to  support  our  claims 
about  the  advantages  of  such  a  process.  As  a  separate  issue  we  will  address  the  unique  problem 
of  the  Block  II/IIA  eclipsing  satellites  and  develop  a  special  solar  pressure  model  for  them. 
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Orbital  geometry  and  coordinate  systems 

The  GPS  satellites  occupy  six  nearly-circular  orbital  planes  with  radius  of  approx.  26,500  km 
and  inclination  of  55°.  The  attitude  of  the  satellites  is  determined  by  the  requirement  that  the 
navigation  antennae  along  the  spacecraft-fixed  Z  axis  point  toward  the  geocenter  and  that  the 
solar  array,  mounted  along  the  spacecraft-fixed  Y  axis,  be  facing  the  Sun.  These  two 
requirements  force  the  spacecraft  to  yaw  constantly,  keeping  the  spacecraft-fixed  X  and  Z  axes 
always  in  the  plane  defined  by  the  Earth,  spacecraft  and  Sun.  The  X  axis  is  required  to  lie  in  the 
half  plane  containing  the  Sun  and  this,  except  at  the  two  singular  points  when  the  Earth,  Sun  and 
spacecraft  are  co-linear,  defines  the  attitude  of  the  spacecraft  uniquely.  If  the  spacecraft  is 
assumed  to  be  symmetric  around  the  XZ  plane  (which  is  the  assumption  built  into  the  ROCK 
models  (Porter,  1983))  then  the  solar  pressure  is  a  two-dimensional  vector  function  of  a  single 
variable,  s,  the  angle  between  the  spacecraft-Sun  direction  and  the  spacecraft-fixed  Z  axis 
(Figure  1).  By  definition,  s  is  allowed  to  vary  in  [0°,180°],  but  its  range  during  a  full  orbit 
revolution  depends  on  the  value  of  the  orbit's  beta  angle.  The  beta  angle  is  the  angle  between  the 
spacecraft-Sun  direction  and  the  orbital  plane  (Figure  2).  It  is  defined  to  be  positive  if  the 
spacecraft-Sun  direction  forms  an  acute  angle  with  the  orbit  normal  and  negative  otherwise.  By 
definition,  beta  ranges  in  [-90°,90°]  but  due  to  the  particular  inclination  of  the  GPS  orbits  it  is 
limited  to,  approximately,  [-75°, 75°].  The  range  of  e  at  a  given  beta  is  [1(31,  180°-I(3I],  It  follows 
that  for  high  beta  angles  the  solar  pressure  hardly  varies  while  for  low  beta  angles  the  solar 
pressure  goes  through  more  of  its  possible  range. 


SUN  w 


Figure  1.  The  XYZ  and  UVW  coordinate  systems.  The  orbital  geometry  corresponds  to  (3  =  0°  as 
viewed  from  above  the  orbital  plane. 
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Figure  2.  The  orbital  geometry  and  the  XYZ  coordinate  system. 


The  spacecraft-fixed  XYZ  coordinate  system  is  useful  for  expressing  the  solar  pressure  and, 
indeed,  the  ROCK  and  the  T20  models  are  expressed  in  this  coordinate  system.  Another  useful 
coordinate  system  is  the  UVW  system  defined  as  follows:  U  is  pointing  along  the  Sun-spacecraft 
direction,  V  is  pointing  along  the  cross  product  between  the  U  direction  and  the  spacecraft- 
geocenter  direction  (and  is,  therefore,  identical  to  the  Y  coordinate)  and  W  completes  the  system 
(Figure  1).  Since  the  bulk  effect  of  the  solar  pressure  is  to  "push"  the  spacecraft  away  from  the 
Sun,  it  is  expected  that  the  dominant  component  of  the  pressure  will  be  in  the  U  direction. 
Consequently,  in  this  system  the  solar  radiation  attains  a  form  which  is  closest  to  a  one¬ 
dimensional  function  and  it  is,  therefore,  most  natural  to  express  the  solar  radiation  in  this 
system. 


The  approach  for  non-eclipsing  satellites 

We  approach  eclipsing  and  non-eclipsing  satellites  differently.  For  the  non-eclipsing  satellites  we 
use  a  fully  empirical  model,  the  parameters  of  which  are  estimated  based  on  a  least  square  fit  to 
“truth”  orbits.  For  eclipsing  satellites  we  used  a  hybrid  approach  where  a  new  model,  a  variant  of 
the  T20  model,  termed  here  the  “split  T20”  model,  was  derived  from  first  principles  and  a 
limited  set  of  model  parameters  was  tuned  manually.  Here  are  the  main  elements  of  our 
approach  for  non-eclipsing  satellites. 

We  express  the  solar  radiation  model  as  a  truncated  Fourier  expansion.  The  expansion  variable 
was  chosen  to  be  £  -  the  Earth-satellite-Sun  angle  -  since,  for  GPS  satellites  following  their 
nominal  attitude  algorithm,  the  satellite-Sun  geometry  is  only  a  function  £.  The  expansion  is  fully 
parameterized,  i.e.,  the  expansion  coefficients  are  presumed  unknown  a-priori.  A  satellite 
trajectory  is  then  created,  based  on  the  empirical  solar  radiation  model,  and  the  model  parameters 
are  adjusted  iteratively  to  fit  the  trajectory  to  a  multi-day  “truth”  orbit  in  a  least  square  sense.  The 
long-arc  “truth”  orbit  is  a  key  element  in  our  approach.  It  is  the  concatenation  of  multiple  daily 
“truth”  orbits  and  it  allows  us  to  retrieve  information  about  the  satellite’s  dynamics  that  is  not 
present  in  the  individual  daily  “truth”  orbits.  The  following  example  illustrates  this  crucial 
concept:  Imagine  having  to  estimate  the  mass  of  a  point-like  planet  around  which  a  satellite  is  in 
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an  8-day-period  circular  orbit.  The  input  is  a  set  of  one-day-long  “truth”  ephemerides  that  were 
estimated  based  on  some  tracking  information  and  the  (erroneous)  assumption  that  the  satellite  is 
free  of  forces  (i.e.,  the  satellite  has  constant  velocity).  As  a  result  of  mismodeling  the  dynamics, 
each  one-day  ephemeris  is  a  straight  line,  but  because  they  are  being  constrained  by  some 
tracking  information  they  cannot  deviate  significantly  from  the  actual  circular  orbit. 
Concatenating  them  will  create  a  long-arc  “truth”  orbit  which  resembles  a  circle  and  allows  for 
the  resolution  of  the  parameters  of  a  proper  model  of  the  dynamics  (see  Figure  3).  By  this 
analogy,  solar  radiation  mismodeling  errors  that  were  built  into  the  daily  GPS  “truth”  solutions 
can  be  backed  out  if  we  concatenate  a  sufficient  number  of  these  daily  orbits  and  use  a  properly 
parameterized  model  to  represent  the  solar  radiation  forces. 


Figure  3.  Stringing  together  straight  segments  reveals  the  underlying  circle. 


In  our  Block  II/IIA  prototype  system  we  chose  the  daily  GPS  orbit  solutions  that  are  produced  at 
JPL  for  the  IGS  as  our  short-arc  “truth”  orbits.  With  a  3-dimensional  error  of  about  15  cm  RMS, 
these  orbits  are  of  the  highest  accuracy  currently  available  (Watkins  et  al.,  1997).  They  are  based 
on  the  T20  solar  pressure  model  but  employ  the  JPL-developed  reduced  dynamics  technique  to 
mitigate  the  short-comings  of  the  T20  model  and  other  models  by  estimating  small  stochastic 
accelerations  (Zumberge  et  al.,  1995). 

Observation  of  the  GPS  navigation  signal  are  made  in  an  Earth-fixed  system  and  the  most 
accurate  representation  of  the  orbits  is  available  in  that  system.  Since  orbit  integration  is  carried 
out  most  conveniently  in  an  inertial  frame  we  have  to  transform  the  solution  from  Earth-fixed 
coordinates  to  inertial  coordinates.  To  minimize  the  aliasing  of  Earth-orientation  errors  into  solar 
pressure  parameters  we  used  the  VLBI-based  IERS  Bulletin-B  values  of  Earth-orientation 
parameters  to  transform  the  solution  between  the  two  frames. 

The  daily  JPL  solutions  span  30  hours,  centered  around  noon  of  each  day.  To  form  a  set  of  daily 
solutions  that  can  be  smoothly  concatenated  we  cosine-average  (see  Appendix  A)  each  pair  of 
consecutive  daily  ephemerides  and  extracted  the  central  24  hours.  We  eliminated  all  satellites 
that  are  known  to  be  problematic  (SVNs  16,  23,  30,  33,  40)  and  those  satellite-days  containing 
maneuvers,  for  the  cleanest  possible  data  set.  The  data  span  is  from  July,  1995,  to  May,  1996. 
The  set  of  transformed,  smoothed,  truncated  and  cleaned  daily  ephemerides  form  our  “truth” 
orbits.  In  the  non-eclipsing  case  we  also  eliminated  all  eclipsing  satellites  from  the  data  set. 
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We  can  now  form  very  long  multi-day  arcs  by  concatenating  consecutive  daily  ephemerides.  Our 
approach  calls  for  estimation  of  model  parameters  based  on  a  least-square  fit  to  a  multi-day 
“truth”  ephemeris.  It  is  clear  that  the  longer  the  arc  the  more  information  about  the  dynamics  can 
be  retrieved  from  it  and  many  months  of  data  are  needed  the  resolve  the  fine  details  of  the  solar 
radiation  model.  Unfortunately,  it  is  not  practical  to  integrate  GPS  orbits  more  than  a  few  days 
because  round-off  error  will  tend  to  contaminate  the  integration  of  the  state  and  variational 
partials.  We  found  it  most  practical  to  perform  10-day  integrations.  Hence,  we  had  to  invent  a 
method  to  combine  estimates  based  on  many  10-day  fits.  This  is  not  a  trivial  problem  because  of 
the  non-linear  nature  of  the  estimation  problem.  This  non-linearity  also  makes  it  difficult  to 
combine  estimates  of  model  parameters  for  different  satellites  into  one. 


The  non-linear  estimation  problem 

Nominally,  we  assume  that  all  Block  II  and  Block  IIA  satellites  should  have  identical  solar 
pressure  models.  This  stems  from  the  assumption  that  they  all  have  the  same  structure  and  obey 
the  same  attitude  laws.  But  the  actual  solar  pressure  is  a  function  of  the  mass  of  the  satellite  and 
of  the  solar  flux  constant.  The  mass  of  each  satellite  is  different  and  both  the  mass  and  the  solar 
flux  constantly  change  in  time.  This  makes  it  difficult  to  combine  estimates  across  different 
satellites  and  different  data  arcs.  Such  a  combination  poses  a  non-linear  estimation  problem.  For 
simplicity,  assume  that  the  solar  radiation  force  is  one-dimensional  and  is  given  by: 

F  =  X  (B0  +  B,  cos  s  +  A,  sin  e  +B2  cos  2e  +  A2  sin  2s  +  ...) 

where  X  is  a  function  of  the  satellite’s  mass  and  of  the  solar  flux  constant.  Estimating  X,  as  well 
as  B0,  B , ,  etc.,  from  a  fit  to  the  10-day  arcs  does  not  always  converge,  partially  because  the 
expansion  coefficients  are  considered  completely  unknown  a-priori.  Also,  there  is  no  unique 
solution  to  this  problem.  Instead,  we  chose  to  estimate  X  B0,  X  B,,  etc.,  thus  postponing  the  non¬ 
linear  problem  to  the  combination  phase.  In  this  phase  we  need  to  combine  the  estimates  X,  B0,  X; 
B,,  ...  from  different  satellites  and  different  time  periods  into  a  more  robust  estimate  of  B0,  B,, 
etc. 

To  describe  our  solution  to  this  non-linear  estimation  problem  we  need  a  more  general  notation. 
Let  a“  be  a  set  of  estimated  parameters  where  i  =  l,...,n  labels  parameters  of  the  solar  pressure 
model  and  a  =  l,...,m  labels  distinct  10-day  batches  and  satellites. 

We  also  have  the  covariances  (using  summation  notation): 

<a“,a/>  =  6apCaij. 

We  associate  a  scale  factor,  Xa,  with  each  set  of  estimated  parameters,  expecting  that 

a“=  Xa  a;. 

We  then  estimate  Xa  and  a;  to  minimize 

J  =  1/22(a“-Xaai)(Ca1)ij(aJ“-XaaJ), 

where  the  explicit  sum  is  over  a.  Again,  this  is  a  non-linear  problem  but  notice  that  Xa  can  be 
solved  analytically  as  a  function  of  a;  from  the  m  equations:  dJ/dXa  =  0.  When  this  is  substituted 
back  into  dJ/da;  =  0  we  obtain  a  system  of  n  (highly)  non-linear  equations  for  the  n  unknowns  a;. 
We  solve  this  system  numerically.  The  redundancy  of  the  system  is  eliminated  by  fixing 

a,  =  H|a;"|. 

With  this  choice  our  numerical  scheme  always  converged.  It  is  not  important  that  the  estimated 
values  of  a;  obtained  in  this  scheme  contain  an  arbitrary  scale  since  any  precise  realization  of  a 
GPS  solar  radiation  model  must  include  the  estimation  of  a  scale  factor.  The  only  thing  that  is 
required  of  the  estimates,  a;,  is  to  be  “reasonable”  and  the  choice  of  a,  =  E[a“]  seems  to 
guarantee  that.  As  a  sanity  check  we  verified  that  the  estimates  for  Xa  are  all  within  a  few  percent 
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of  1,  in  agreement  with  our  knowledge  of  time- variations  in  the  solar  flux  constant  and  variations 
in  mass  properties  between  satellites.  There  is  no  practical  limit  on  the  number  of  satellites  or  10- 
day  batches  that  can  be  combined  with  this  scheme. 

Evaluation  of  solution  quality  is  the  last  step  in  our  approach.  This  is  done  by  testing  each 
model’s  ability  to  best-fit  a  pre-selected  set  of  truth  ephemerides.  There  are  four  types  of  tests: 

1.  best  fit  to  single-day  “truth”  arcs, 

2.  best  fit  to  4-day  “truth”  arcs, 

3.  best  fit  to  10-day  “truth”  arcs, 

4.  fit  to  “truth”  during  a  4-day  prediction  period  based  on  best-fit  to  previous  4-day  “truth” 
orbit. 

Each  type  of  test  consists  of  many  (MOO)  satellite-arcs  to  provide  a  robust  evaluation  of  the 
model’s  performance.  Although  the  model’s  parameters  are  fixed  to  their  estimated  value  during 
the  test  run  there  is  a  set  of  auxiliary  parameters  that  are  always  being  estimated  to  provide  best 
fit,  in  addition  to  the  six  epoch  state  parameters.  Not  least  among  these  is  the  overall  scale  factor 
which,  as  mentioned  above,  always  needs  to  be  estimated  to  account  for  the  small  but 
undetermined  scale  that  is  present  in  the  model’s  estimated  parameters.  Other  auxiliary 
parameters  vary  depending  on  the  type  of  the  model  being  tested  but  they  usually  include  the  Y 
bias  and,  perhaps,  one  or  two  additional  parameters  that  are  deemed  not  constant  in  time  or 
between  satellites.  The  significance  of  the  various  test  types  increases  in  the  order  they  are  listed 
above.  Since  we  usually  adjust  at  least  8  auxiliary  parameters  (six  epoch  state,  overall  scale  and 
Y  bias)  our  ability  to  fit  short  orbit  arcs  is  largely  affected  by  errors  in  our  “truth”  ephemerides 
rather  than  the  quality  of  the  estimated  model.  The  test  that  is  most  sensitive  to  the  model’s 
quality  is  the  prediction  fit.  In  this  fit  we  fix  the  value  of  the  auxiliary  parameters  to  those 
estimated  for  best  fit  in  the  4-day  test  case.  We  then  continue  to  integrate  the  orbit  for  four  more 
days  and  compare  the  last  4  days  of  the  trajectory  to  the  “truth”.  The  intention  here  is  not  to 
generate  the  best  prediction  possible  (this  will  require  a  different  estimation  scheme  during  the 
first  four  days)  but,  rather,  to  compare  the  performance  of  various  models  in  a  relative  sense. 

In  summary,  our  approach  consists  of  4  steps: 

1.  obtaining  daily  “truth”  orbits  spanning  many  days  and  forming  10-day  arcs  by 
concatenation, 

2.  for  each  10-day  arc  and  satellite  adjusting  the  parameters  of  a  Fourier-like  solar  radiation 
model  to  best  fit  the  “truth”  orbit, 

3.  combining  the  estimates  from  each  10-day  arc  and  for  each  satellite,  together  with  their  full 
covariance  matrix  into  a  single,  robust  estimate  for  the  model’s  parameters, 

4.  evaluating  the  resulting  model  (with  the  estimated  parameters)  by  comparing  against 
“truth”  ephemerides. 


Results  for  non-eclipsing  satellites 

The  model  space  for  the  possible  solar  radiation  model  is  infinite.  For  the  purpose  of  this 
prototype  development  we  have  limited  our  search  to  very  low  order  Fourier  expansions.  We 
chose  the  UVW  coordinate  system  for  the  expansion  because  we  expect  the  fewest  Fourier  terms 
in  this  system.  We  postulate  the  expansion: 

Fu  =  f(m,r)  s  (C_U0  +  C_U,  cos  8  +  S_U,  sin  8  +  C_U2  cos  28  +  S_U2  sin  28  +  ...) 

Fv  =  C_V0  +  C_Y,  cos  8  +  S_V,  sin  8  +  ... 

Fw  =  f(m,r)  s  (C_W,  cos  8  +  S_W,  sin  8  +  C_W2  cos  2e  +  S_W2  sin  28  +  ...) 
where  f(m,r)  is  a  function  of  the  satellite’s  mass  and  distance  from  the  Sun.  s  is  an  overall  scale 
factor  that  is  set  to  1  during  the  estimation  of  the  model  parameters.  Notice  that  the  V  structure  is 
different  than  the  U  and  W  structures.  That  is  because  the  nominal  symmetry  of  the  spacecraft 
about  the  UW  plane  leads  to  consideration  of  forces  orthogonal  to  this  plane  as  anomalous  and 
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fewer  terms  are  expected  in  the  V  structure.  Indeed,  until  this  study  revealed  energy  in  higher 
harmonics,  C_V0  was  considered  the  only  significant  component  (since  the  V  coordinate  is 
identical  to  the  Y  coordinate,  C_V0  is  just  another  name  for  the  Y  bias).  We  experimented  with  a 
few  low  order  truncations  of  the  above  expansion.  For  a  given  truncation,  all  of  the  harmonics 
are  estimated  using  10-day  batches  spanning  9  months  and  all  “good”,  non-eclipsing  satellites. 
Since  the  overall  scale,  s,  was  set  to  1,  these  estimates  effectively  contain  an  unknown  scale.  The 
scale  may  vary  between  different  satellites  and  different  10-day  arcs  but  is  the  same  for  all 
estimated  harmonics  in  a  given  10-day  arc  and  a  given  satellite. 

Figures  4-6  illustrate  the  range  of  the  estimated  values,  together  with  their  formal  error  for  a  few 
of  the  most  interesting  harmonics.  These  figures  demonstrate  the  subtleties  of  the  estimation 
process.  For  low  beta  angles,  the  U  coordinate  lies  close  to  the  orbital  plane  and  the  Fu 
component  of  the  solar  pressure  is  dynamically  significant  (Colombo,  1992).  As  the  beta  angle 
increases,  the  dynamic  significance  of  the  Fu  harmonics  diminishes.  This  is  reflected  by  the  high 
formal  errors  for  C_U0  as  beta  increases  in  absolute  value  (Figure  4).  Figure  5  depicts  the 
estimated  values  of  C_V0.  For  high  beta  angles  the  V  coordinate  points  along-track  and  is, 
therefore,  dynamically  very  significant.  For  low  beta  angles,  the  V  coordinate  rotates  about  the 
radial  direction  (the  Z  coordinate),  and  because  it  does  not  point  strictly  across-track  it  is  also 
dynamically  significant.  This  fact  is  reflected  as  low  formal  errors  for  all  V  harmonics  at  all  beta 
values.  One  of  the  most  interesting  results  of  this  study  is  the  discovery  of  significant  energy  in 
the  first  cosine  harmonic  in  the  V  direction.  Figure  6  depicts  the  estimated  values  of  this 
harmonic.  It  shows  a  remarkable  coherence  in  the  estimated  value  between  all  satellites  and  a 
very  low  level  of  formal  error.  Plots  of  other  parameters  are  presented  in  appendix  B. 
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Figure  4.  Estimates  of  C_U0,  with  formal  errors  for  non-eclipsing  (1(11  >  14)  Block  IIA  satellites. 
The  estimates  for  high  beta  angles  are  usually  off  the  scale  of  this  plot. 


•  GPS22  •  GPS25  •  GPS27  •  GPS29  •  GPS32  •  GPS35  •  GPS37 

•  GPS24  ®  GPS26  •  GPS28  •  GPS31  •  GPS34  •  GPS36  ®  GPS39 


Figure  5.  Estimates  of  C_V0  (Y  bias),  with  formal  errors  for  non-eclipsing  (1(11  >  14)  Block  IIA 
satellites  .  The  formal  errors  are  usually  too  small  to  be  seen. 
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Figure  6.  Estimates  of  C_V\,  with  formal  errors  for  non-eclipsing  (1(31  >  14)  Block  IIA  satellites  . 
The  formal  errors  are  usually  too  small  to  be  seen. 

After  we  obtain  the  set  of  estimated  values  for  each  harmonic  from  all  satellites  and  10-day  arcs, 
we  can  combine  them  to  form  a  single,  robust  estimate  for  each  harmonic.  At  this  point  we  need 
to  select  the  set  of  harmonics  that,  first,  are  assumed  constant  after  removing  the  satellite-  and 
arc-specific  scale  and,  second,  are  dynamically  important.  Here  some  trial  and  error  is  required 
to  pick  a  good  set  of  harmonics  to  determine  our  model.  We  have  chosen  not  to  merge  the  C_V0 
estimates  because  they  are  traditionally  always  being  estimated  and,  being  anomalous,  we 
suspect  that  they  might  not  be  scaled  the  same  way  as  the  U  and  W  harmonics.  The  estimates  of 
C_V,  clearly  demonstrated  that  this  harmonic  is  not  a  constant  and  it  was  also  omitted  from  the 
combination  process.  After  some  experimentation  the  following  set  of  harmonics  was  chosen  for 
the  combination  process: 

C_U0,  C_Uh  S_U,,  S_U2,  S_W , ,  C_W,  and  S_W2.  Although  not  constant,  our  analysis  showed 
that  the  presence  of  the  C_Vi  harmonic  improves  significantly  the  quality  of  the  test  fits. 
Consequently,  we  approximated  the  estimated  values  of  C_V,  as  the  following  function  of  beta 
(see  Figure  7): 

C_V,  =  0.1  +  0.5  sin  (3  +  0.3/sin  (3 

The  impact  that  the  presence  of  C_V,  makes  on  the  quality  of  the  model  is  illustrated  in  Figure  8. 

Our  final  model  takes  the  form  (full  precision  and  formal  errors  are  reported  in  appendix  C): 

Fu  =  f(m,r)  s  (8.98  -  0.15  cos  e  +  0.71  sin  e  +  0.1  sin  2s) 

Fv  =  v_bias  +  (0.1  +  0.5  sin  (3  +  0.3/sin  (3)  cos  s 
Fw  =  f(m,r)  s  (-0.03  sin  e  +  0.73  cos  s  -  0.74  sin  2s) 

were  s  and  v_bias  are  always  estimated  following  the  conventional  practice.  This  model  is 
labeled  GSPM.II.97+CY,.  If  the  C_V,  harmonic  is  excluded  the  resulting  model  is  labeled 
GSPM.II.97.  This  model  is  as  simple  in  structure  as  the  T20  model  and  it  has  only  a  few  more 
coefficients.  The  T20  model  is: 

Fx  =  f(m,r)  s  (-8.96  sin  e  +  0.16  sin  3e  +  0.1  sin  5e  -  0.07  sin  7e) 

Fy  =  y_bias 

Fz  =  f(m,r)  s  (-8.43  cos  e) 


Figure  7.  Fitting  a  functional  form  to  C_V,. 
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Figure  8.  Effects  of  the  C_V,  harmonic  on  the  ability  of  the  model  to  fit  the  1-4-  and  10-day 
“truth”  ephemerides.  Top:  the  GSPM.  11.97  model  with  additional  estimation  of  solar  scale  and 
C_V0  (Y  bias).  Bottom:  same  as  in  top  panel  but  with  the  addition  of  C_V,  to  the  model. 


The  performance  of  our  new  model  is  summarized  in  Figure  9.  Most  of  the  improvement  in 
model  performance  is  due  to  the  presence  of  the  C_V,  harmonic.  But  even  if  this  harmonic  is 
omitted,  our  flight-data-based  model  performs  as  well  as  or  better  than  the  T20  model. 

We  have  experimented  with  various  other  low  order  Fourier  expansions  expressed  both  in  the 
UVW  and  XYZ  coordinate  systems.  The  differences  in  overall  accuracy  were  found  to  be 
negligible  and,  consequently,  they  are  not  presented  here. 
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10  days  4  day  prediction 


Figure  9.  Performance  of  the  flight-data-based  model  in  comparison  to  the  T20  model.  The 
performance  is  measured  in  terms  of  mean  RMS  fit  to  “truth”  ephemerides.  In  all  cases  solar 
pressure  and  Y  bias  (C_V0)  are  the  only  estimated  parameters. 
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Approach  for  eclipsing  satellites 

The  problem  in  modeling  solar  pressure  for  eclipsing  GPS  satellites  stems  from  the  fact  that  GPS 
Block  II/IIA  satellites  display  a  non-nominal  attitude  behavior  during  eclipse  season.  The  non- 
nominal  behavior  is  limited  to  the  yaw  axis  and  to  regions  of  the  orbit  around  “noon”  and 
“midnight”.  Orbit  noon  is  defined  as  the  point  on  the  orbit  that  is  closest  to  the  Sun.  Orbit 
midnight  is  defined  as  the  point  on  the  orbit  that  is  farthest  from  the  sun.  The  nominal  attitude 
calls  for  the  satellite’s  Z  axis  (the  yaw  axis)  to  point  toward  the  geocenter  and  for  the  satellite  X 
axis  to  lie  in  the  Sun-Earth-satellite  plane  and  point  in  the  general  direction  of  the  Sun.  Also,  the 
large  solar  array  is  assumed  to  be  directly  facing  the  sun.  The  T20  model  erroneously  assumes 
that  GPS  satellites  maintain  nominal  attitude  at  all  times.  The  non-nominal  behavior  manifests  as 
a  non-nominal  pointing  of  the  X  axis  which,  in  turn,  leads  to  mispointing  of  the  solar  array.  This 
happens  during  every  shadow  crossing  and  up  to  45  minutes  thereafter,  and  in  the  vicinity  of 
orbit  noon  at  low  beta  angles  (<  5°)  for  up  to  30  minutes.  This  non-nominal  pointing  is  the 
largest  error  source  in  the  solar  pressure  model  for  eclipsing  satellites.  (For  background  on  the 
eclipsing  problem  please  refer  to  Bar-Sever  et  al.,  1996  and  Bar-Sever,  1996).  The  actual  attitude 
of  an  eclipsing  satellite  can  be  modeled  during  all  phases  of  the  orbit  and  is  described  in  Bar- 
Sever  (1996).  The  actual  attitude  model  has  been  part  of  the  GIPSY-OASIS  software  since 
September,  1994,  but  the  lack  of  a  compatible  solar  pressure  model  prevented  its  use  in  the 
dynamic  modeling  of  the  satellite’s  trajectory  until  now. 

Our  approach  for  developing  a  solar  pressure  model  for  Block  II/IIA  eclipsing  satellites  differs 
from  our  approach  for  non-eclipsing  satellite  for  the  following  reasons: 

1.  the  orbital  geometry  of  eclipsing  satellites  is  non-periodic, 

2.  future  GPS  blocks  (Block  HR  and  Block  IIF  satellites)  will  behave  differently  during 
eclipse  periods, 

3.  the  errors  in  the  current  model  (T20)  are  much  larger  for  eclipsing  satellites  than  they  are 
for  non-eclipsing  satellites. 

Item  number  1  suggests  that  the  Fourier  expansion  approach  would  be  much  more  difficult  to 
implement  for  eclipsing  satellites,  requiring  possibly  another  expansion  variable.  Item  number  2 
suggests  that  methods  used  in  developing  solar  pressure  for  Block  II/IIA  eclipsing  satellites 
cannot  be  readily  used  for  Block  HR  and  IIF  eclipsing  satellites.  Item  number  3  suggests  that 
even  a  crude  model  might  improve  upon  the  currently-available  model.  Consequently,  we 
approached  the  problem  of  improving  the  solar  pressure  model  for  eclipsing  satellites  not  as  a 
prototype  for  future  implementation  but  as  a  means  for  itself.  To  this  end  we  proceeded  to 
develop  a  force  model  from  first  principles. 

We  assume  that  when  the  satellite  is  nominally  oriented  we  can  use  the  nominal  solar  pressure 
model.  This  can  be  the  T20  model  or  one  of  our  newly  developed  models.  For  simplicity  we 
based  this  development  on  the  T20  model.  We  further  assume  that  during  shadow  crossing  no 
radiation  forces  are  acting  on  the  satellite.  What  is  left,  is  to  design  a  model  which  will  be  valid 
during  the  non-nominal  pointing  periods,  which  occur  mainly  after  shadow  exit.  During  this 
post-shadow  maneuver  the  satellite  is  yawing  rapidly  to  regain  its  nominal  yaw  attitude  which 
was  lost  during  shadow  crossing  (Bar-Sever,  1996). 

We  termed  the  new  model  “split  T20”  for  reasons  that  will  be  obvious  soon.  To  derive  it  we 
assumed  that  the  total  force  indicated  by  the  T20  model  is  a  superposition  of  two  contributions:  a 
radiation  force  contributed  by  the  solar  array  and  a  radiation  force  contributed  by  the  main  body 
of  the  spacecraft.  To  determine  the  two  contributions  we  proceeded  as  follows.  Let  Fx  and  Fz  be 
the  X  and  Z  components  of  the  T20  force  in  the  X  and  Z  coordinates.  In  the  UVW  coordinate 
system  the 
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T20  force  has  components  Fu  and  Fw  were, 

Fu  =  -Fx  sin  8  -  Fz  cos  8  =  8.43  +  f(e) 

Fw  =  Fx  cos  8  -  Fz  sin  8  =  g(e) 

We  focus  on  the  constant  contribution  to  Fu.  Solar  array  radiation  forces  are  nominally  constant 
because  the  solar  array  is  nominally  always  facing  the  Sun.  Hence  the  solar  array  radiation  force 
constitutes  part  of  the  8.43  constant.  But  the  body  radiation  force  should  also  have  a  zero-order 
term  that  will  contribute  toward  the  8.43  constant.  We  now  introduce  a  new  parameter,  r,  equal 
to  the  fraction  the  solar  array  contribution  is  from  8.43.  Based  on  available  information  about  the 
satellite  structure,  the  value  of  r  should  be  approximately  0.9.  Given  r,  we  have  split  the  T20 
model  into  its  two  components:  a  constant  solar  array  radiation, 

Fu  =  r8.43, 

and  the  remainder,  due  to  main  body  radiation  ,  given  by 
Fu  =  (1  -  r)  8.43  +  f(e) 

Fw  =  g(e). 

Once  we  split  the  T20  model  into  its  two  components  we  make  the  following  additional 
assumptions:  1.  The  satellite’s  body  is  XY- symmetric  and,  2.  The  radiation  force  on  the  solar 
panels  obeys  the  rule: 

F  =  -C  cos  a  [(  2/3  6  +  2  p  cos  a)  N  +  (1  -  p)  S)] 

where  C  is  a  function  of  the  area  to  mass  ratio,  a  is  angle  between  the  spacecraft-Sun  direction,  S 
and  N  -  the  normal  to  the  solar  array.  6  and  p  are  the  diffusive  and  specular  reflectivity, 
respectively,  of  the  solar  array.  Given  these  assumption  we  simply  rotate  the  main-body-induced 
radiation  force  from  its  nominal  direction  around  the  Z  axis  by  the  amount  of  yaw  error.  The 
solar  array  contribution  is  now  split  along  the  S  and  N  directions  according  to  the  above  formula. 
C,  which  is  not  known  a-priori  is  determined  by  setting  N  =  S  and  F  =  r  8.43  S.  Notice  that 
during  nominal  orientation  N  =  S  and  the  “split  T20”  model  reduces  to  the  T20  model. 

The  biggest  challenge  here  is  the  proper  modeling  of  N,  the  solar  panel  normal,  as  a  function  of 
time.  Following  the  principles  of  the  post-shadow  maneuver  outlined  in  Bar-Sever  (1996),  we 
added  a  model  to  our  GIPSY  software  that  accounts  for  the  motion  of  the  solar  array,  assuming  a 
pitch  rate  of  0.25  deg/sec.  A  critical  element  in  this  model  is  the  determination  of  the  sense  of 
yaw  during  the  post-shadow  maneuver.  As  discussed  in  Bar-Sever  (1996),  this  can  be  very 
difficult  whenever  the  yaw  error  after  shadow  crossing  approaches  180°.  This  issue  has  not  been 
resolved  yet  and,  consequently,  the  performance  of  the  model  is  sub-optimal.  We  estimate  that 
about  10%  of  the  post-shadow  maneuvers  are  modeled  incorrectly  as  a  result  of  this  deficiency. 
This  is  still  a  significant  improvement  over  the  100%  of  post-shadow  maneuvers  mismodeled  by 
the  T20  model. 


Results  for  eclipsing  satellites 

The  model  has  three  parameters:  r,  6  and  p.  Their  values  were  determined  manually  using  trial 
and  error  so  as  to  achieve  best  fit  to  the  “truth”  ephemerides.  The  chosen  values  are: 
r  =  0.9. 
p  =  0.25. 

6  =  0.3. 

The  performance  of  the  “split  T20”  model  was  compared  to  that  of  the  T20  model  in  the  same 
way  as  was  done  for  the  non-eclipsing  satellites,  namely,  by  comparing  their  ability  to  fit  “truth” 
ephemerides  of  various  lengths.  The  results  as  summarized  in  Figure  10,  demonstrate  clearly  the 
very  significant  improvement  realized  by  the  “split  T20”  model.  In  this  comparison  C_Y,  (C_V,) 
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is  being  estimated  in  all  cases  in  addition  to  the  solar  scale  and  Y  bias.  The  singular  behavior 
exhibited  by  C_Y,  at  low  beta  angles  (Figure  7)  precludes  us  from  using  its  functional  form  as  a 
fixed  part  of  the  model,  as  we  did  for  the  non-eclipsing  satellites. 


1  day 


4  days 


Figure  11.  Performance  of  the  “split  T20”  model  in  comparison  to  the  T20  model.  The 
performance  is  measured  in  terms  of  mean  RMS  fit  to  “truth”  ephemerides.  In  all  cases,  overall 
scale,  Y  bias  (C_Y0)  and  C_Y,  are  the  only  estimated  parameters. 
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Discussion 

We  have  demonstrated  here  the  successful  development  of  a  flight-data-based  solar  pressure 
model  for  non-eclipsing  GPS  Block  II/IIA  satellites.  Although  the  scope  of  this  work  did  not 
allow  for  an  exhaustive  search  for  the  best  possible  model,  a  significant  improvement  over  the 
T20  model  was  achieved  with  the  limited  set  of  models  we  tested.  The  models  are  all  in  the  form 
of  a  Fourier  expansion  in  8.  Almost  all  of  the  improvement  in  model  quality  comes  from  a 
newly-found  structure  in  the  Y  (V)  direction  which  we  were  able  to  express  as  a  simple  function 
of  8  and  (3.  This  structure  can  be  easily  embedded  in  any  existing  solar  pressure  model  to  yield 
significant  improvement.  The  term  “solar  radiation”  should  not  be  construed  as  accounting  only 
for  radiation  emanating  directly  from  the  Sun.  In  fact,  the  flight-data-based  model  cannot 
distinguish  between  Earth  radiation  in  the  optical  (albedo)  and  infra-red  range  and  direct  solar 
radiation.  Instead,  it  lumps  them  together,  which  is  part  of  its  strength.  The  term  “solar  radiation 
model”  is  still  appropriate  since  the  ultimate  source  for  all  radiation  effects  is  the  Sun. 

In  terms  of  overall  accuracy  we  have  not  found  any  significant  difference  between  the  two 
coordinate  systems  we  experimented  with,  although,  as  expected,  the  UVW  system  was  more 
“economical”,  requiring  fewer  terms  in  the  expansion.  It  is  noteworthy,  though,  that  given  our 
limited  model  space  we  found  it  difficult  to  improve  significantly  upon  the  T20  model’s  XZ 
components  but,  also  noteworthy,  is  the  ease  with  which  we  could  get  a  similar  level  of 
accuracy. 

It  is  evident  from  the  level  of  error  still  present  in  orbit  prediction  (Figure  9)  that  there  is  room 
for  further  modeling  improvements.  Since  our  model  space  in  this  study  had  to  be  limited  to  low 
order  Fourier  expansions  we  might  conclude  that  further  improvements  can  be  realized  by 
incorporating  higher  order  terms.  The  significant  magnitude  of  the  newly  discovered  first  order 
term  in  the  Y  direction  seems  to  suggest  the  possibility  of  higher  order  structure  in  that  direction. 
Feeding  more  data  into  our  low  order  models  had  only  a  marginal  effect.  Improvements  may  also 
be  realized  from  tuning  the  model  for  individual  satellites  although  more  data  may  be  needed  to 
compensate  for  the  increase  in  degrees  of  freedom. 

The  ability  of  the  new  models  to  improve  GPS  orbits  (rather  than  fit  them)  from  raw  observations 
still  needs  to  be  demonstrated.  It  is  expected  that  iterating  on  this  process  will  result  in  further 
improvement  to  the  model  and  to  the  model-based  determined  orbits. 

An  important  by-product  of  this  development  effort  has  been  the  a  insight  into  the  dynamics  of 
GPS  satellites  afforded  by  examining  highly  accurate,  quality-controlled,  long  time  series  of  GPS 
ephemerides  in  a  well  defined  inertial  frame.  To  our  knowledge  this  is  the  first  time  that  such 
high-quality  data  has  been  produced  and  examined.  Our  preliminary  analysis  revealed,  for  the 
first  time,  hard  evidence  relating  the  Y  bias  to  solar  array  misalignment  -  a  long-standing 
hypothesis  (Fliegel  et  al.,  1992).  Figure  5  displays  a  baffling  asymmetry  in  the  behavior  of  the  Y 
bias  estimates  as  a  function  of  beta.  Rearranged  as  a  function  of  time,  a  large  discontinuity  is 
revealed  in  the  estimates  for  seven  satellites,  SVNs  22,  25,  26,  27,32,  35  and  39  (Figure  11).  The 
common  denominator  for  all  these  satellites  is  that  the  yaw  bias  implemented  in  their  attitude 
control  system  changes  sign  simultaneously  with  the  observed  discontinuity.  Prior  to  November, 
1995,  the  GPS  operators  had  routinely  changed  the  sign  of  the  yaw  bias  in  the  satellite  attitude 
control  system  when  the  orbit’s  beta  angle  switched  sign.  Bar-Sever  (1996)  suggested  a  link 
between  this  practice  and  a  simultaneous  jump  that  was  observed  in  the  values  of  the  estimated 
yaw  rates.  From  November  1995  onward,  the  yaw  bias  was  set  permanently  to  +0.5°  and  indeed, 
since  then,  no  discontinuity  was  observed  either  in  the  estimated  yaw  rate  or  in  the  values  of  the 
estimated  Y  bias.  This  connection  suggests  that  the  yaw  bias  is  partially  responsible  for  the 
currently  observed  Y  bias.  It  is  by  no  means  the  only  contributor  to  the  Y  bias  but  it  does  suggest 
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that  even  a  small  amount  of  solar  array  misalignment ,  on  the  order  of  1°,  is  a  likely  cause  of  the 
Y  bias. 


GPS22  — GPS25  — — GPS27  -  GPS29  — • — GPS32  — • — GPS35  •  GPS37 


o  GPS24  — « — GPS26  °  GPS28  °  GPS31  °  GPS34  °  GPS36  -^^GPS39 


Figure  11.  Estimated  values  of  Y  bias  as  a  function  of  time  for  Block  IIA  satellites.  The 
estimates  are  based  on  10-day  fits  to  “truth”.  The  satellites  displaying  discontinuities  across 
eclipse  season  are  denoted  by  the  heavy  lines.  The  satellites  not  displaying  the  discontinuity  are 
denoted  by  the  broken  line.  The  time  span  is  1995  -  1996. 

Significant  improvement  to  the  modeling  of  eclipsing  satellites  was  realized  trough  the  “split 
T20”  model.  This  model  still  uses  the  T20  model  during  the  portion  of  the  orbit  when  the 
satellite  is  oriented  nominally  (about  %85  of  the  time).  Much  of  the  remaining  error  is  due  to 
unresolved  post-shadow  maneuver  ambiguities.  Although  it  is  possible,  in  principle,  to  resolve 
all  the  post-shadow  maneuver  ambiguities  given  enough  flight  data,  it  is  impossible  to  do  so  in 
the  case  of  orbit  prediction,  when  no  such  data  is  available.  This  is  the  main  reason  why  the 
model  for  eclipsing  satellites  will  never  be  as  accurate  as  for  non-eclipsing  satellites. 

One  of  the  perennial  problems  with  GPS  has  been  the  resolution  of  UT1  rate.  In  precise  GPS 
orbit  determination  solar  radiation  parameters  as  well  as  length-of-day  are  being  simultaneously 
adjusted.  Unfortunately,  a  significant  level  of  correlation  between  the  Y  bias  and  UT1  rate 
contributes  to  errors  in  the  realization  of  an  inertial  frame  by  the  GPS.  One  of  the  unique 
advantages  of  our  technology  is  the  ability  to  estimate  Y  bias  in  a  “solid”  inertial  frame  defined 
by  VLBI.  Figure  5  display  these  estimates,  per  satellite,  as  a  function  of  the  beta  angle.  These 
estimates  are  based  on  10-day  fits  over  9  months.  Ignoring  the  discontinuities  explained  earlier, 
the  estimates  exhibit  remarkable  coherence  per  satellite.  These  estimates  can  be  the  basis  of  a  Y 
bias  function  that  is  made  part  of  the  solar  pressure  model.  This  will  eliminate,  or  significantly 
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reduce,  the  need  to  estimate  Y  bias  in  normal  precise  orbit  determination  operations  and  increase 
the  quality  of  the  UT1  rate  estimates. 
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Figure  12.  Estimates  of  model  parameters  based  on  data  arcs  of  various  lengths.  Each  point 
represents  an  estimated  value  of  the  corresponding  parameter  from  a  single  arc.  The  top  panel  is 
for  estimates  based  on  10-day  arcs,  the  middle  panel  is  for  estimates  based  on  30-day  arcs  and 
the  bottom  panel  is  for  estimates  based  on  90-day  arcs. 
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In  order  to  demonstrate  the  usefulness  of  this  technology  for  newly-launched  satellites,  where  the 
amount  of  data  is  very  limited  initially,  we  studied  the  relationship  between  the  quantity  of  the 
data  and  the  quality  of  the  estimates.  To  this  end,  we  set  up  an  experiment  where  we  estimated 
the  model  parameters  based  on  10  days  of  data,  30  days  of  data,  and  90  days  of  data.  The  results 
presented  above  rely,  of  course,  on  9  months  of  data.  In  this  experiment  we  divided  the  9  months 
to  disjoint  sets  of  10-,  30-  and  90-day  arcs  and  we  estimated  the  model  parameters  from  a  fit  to 
each  arc.  The  effect  on  the  estimated  values  is  illustrated  in  Figure  12.  As  can  be  expected,  the 
scatter  in  most  model  parameters  decrease  when  more  data  is  brought  to  bear.  Nevertheless,  if  a 
model  which  is  based  on  a  single  10-day  fit  is  used  during  the  10  day  period  immediately  after 
that,  its  ability  to  fit  a  truth  orbit  is  similar  to  that  of  a  model  based  on  9  months  of  data.  The 
main  advantage  of  having  more  data  is  to  make  the  model  more  robust,  namely,  valid  throughout 
all  mission  phases.  This  result  suggests  that  this  technology  is  suitable  for  use  with  newly 
launched  satellites,  where  the  robustness  of  the  model  can  be  increased  gradually  as  more  data 
become  available,  and  as  preliminary  models  are  refined  in  an  iterative  process. 

An  ideal  opportunity  to  validate  and  to  tune  this  technology  for  newly-launched  satellite  swill 
come  as  Block  HR  satellites  are  about  to  be  launched.  Block  HR  satellites  follow  a  different 
attitude  control  algorithm  than  Block  II/IIA  satellites.  This  requires  some  preparatory  work  to 
model  the  new  algorithm.  We  expect  that  application  of  this  technology  immediately  after  launch 
will  give  rise  to  a  solar  radiation  model  with  much  better  accuracy  than  the  T30  model  (Fliegel 
and  Gallini,  1996)  and,  in  particular,  will  be  able  to  account  properly  for  the  expected  BOL 
anomalies. 

The  ultimate  benefit  from  this  technology  will  be  derived  from  applying  it  to  Block  IIF 
satellites.  Since  the  pre-launch  design  phase  of  the  solar  radiation  model  has  not  begun  yet,  it  is 
possible  to  realize  significant  cost  savings  by  eliminating  it  altogether.  Another  possible  benefit 
for  Block  IIF  satellite  can  be  reaped  if  the  attitude  algorithm  design  process,  currently  taking 
place,  will  avoid  the  limitation  on  accuracy  for  eclipsing  satellites  imposed  by  Block  II,  IIA  and 
IIR  designs.  As  mentioned  above,  the  dynamics  of  Block  II/IIA  (and  we  expect  also  Block  HR) 
eclipsing  satellites  can  never  be  modeled  with  the  same  high  fidelity  as  non-eclipsing  satellites. 
This  is  not  unavoidable,  and  solutions  do  exist  that  eliminate  this  problem. 
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Appendix  A.  Cosine  averaging 

Given  two  30-hour  ephemerides  with  6-hour  overlap,  this  scheme  averages  the  data  points  in  the 
overlap  region  by  giving  them  variable  weight  based  on  their  distance  from  the  edge  of  the 
relevant  30-hour  span.  If  A  is  a  data  point  from  the  last  6  hours  of  the  first  day  and  B  is  the 
corresponding  data  point  from  the  first  6  hours  of  the  second  day,  then  the  cosine  averaging  of  A 
and  B  is 

[A  +  B]  =  A  (1  +  cos(jt  T/6))/2  +B  (1  -  cos(jt  T/6))/2, 
where  T  is  the  time  distance,  in  hours,  between  data  point  A  and  the  beginning  of  the  overlap 
period. 
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Appendix  B.  Estimated  harmonics. 

The  Figures  in  this  Appendix  depict  the  estimated  values  and  formal  errors  for  all  the  parameters 
of  the  GSPM.II.97  model  not  given  above.  The  abscissa  scale  is  always  10-5  Newtons. 
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Appendix  C.  Coefficients  and  formal  errors  of  GSPM.II.97 

C_u0=  8.981090  ±  0.020854 
C_U,=  -0.148668  ±  0.009257 
C_U2  =  -0.000274  ±  0.013199 

SJJ,  =  0.713423  ±  0.032058 

S_U2  =  0.101430  ±  0.005688 

S_Wj  =  -0.025771  ±  0.002933 

S_W2  =  -0.742966  ±  0.010121 

C_W,  =  0.732856  ±  0.014482 

C_W2  =  -0.005984  ±  0.002590 
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